Example: simulation of a satellite in orbit around Earth
The following examples show how to perform a basic orbital simulation: satellite ephemeris, visibility of ground stations, eclipses...
Note: First execute the code below to create the necessary variables (satellite ephemeris, ground stations ...).
// ===================================================== // Satellite ephemeris // The orbit is approximately frozen and Sun-synchronous // Initial mean local time of ascending node (MLTAN) = 10h30 // ===================================================== // Date/time of orbital elements (TREF) t0 = CL_dat_cal2cjd(2013,3,1,12,0,0); // Keplerian mean orbital elements, frame = ECI sma = 7200.e3; // semi major axis ecc = 1.e-3; // eccentricity inc = 98 * %pi/180; // inclination pom = %pi/2; // Argument of perigee mlh = 10.5; // MLTAN (hours) gom = CL_op_locTime(t0, "mlh", mlh, "ra"); // RAAN anm = 0; // Mean anomaly kep0 = [sma; ecc; inc; pom; gom; anm]; // Simulation dates/times (duration = 1 day) t = t0 + (0 : 30/86400 : 1); // Propagate with "lydlp" model (output = osculating elements) kep = CL_ex_propagate("lydlp", "kep", t0, kep0, t, "o"); // Position and velocity in ECI [pos_eci, vel_eci] = CL_oe_kep2car(kep); // Position in ECF pos_ecf = CL_fr_convert("ECI", "ECF", t, pos_eci); // ===================================================== // Other data // ===================================================== // Ground stations geodetic coordinates: // longitude (rad), latitude (rad), altitude (m) sta1 = [1.499*%pi/180; 43.43*%pi/180; 154.0]; sta2 = [-52.64*%pi/180; 5.1*%pi/180; 94.0]; // Earth->Sun and Earth->Moon in ECI Sun_eci = CL_eph_sun(t); Moon_eci = CL_eph_moon(t); |  |  |
1) Check orbit properties
// ===================================================== // Check Sun-synchronicity // ===================================================== // Mean local time of ascending node mlh = CL_op_locTime(t, "ra", kep(5,:), "mlh"); // Plot scf(); plot(t - t0, mlh); xtitle("Mean local time of ascending node (h)", "Days"); CL_g_stdaxes(); | | |
2) Ground station visibility start/end times
// ===================================================== // Ground stations visibility // ===================================================== // Min elevation for visibility min_elev = 10 * %pi/180; // Computation of visibility intervals [tvisi1] = CL_ev_visibilityEph(t, pos_ecf, sta1, min_elev); [tvisi2] = CL_ev_visibilityEph(t, pos_ecf, sta2, min_elev); // Plot scf(); plot2d3(tvisi1(1,:)-t0, (tvisi1(2,:)-tvisi1(1,:))*1440, style=2); plot2d3(tvisi2(1,:)-t0, (tvisi2(2,:)-tvisi2(1,:))*1440, style=5); xtitle("Visibility duration (mn)", "Days"); h = CL_g_select(gca(), "Polyline"); h.thickness = 2; CL_g_stdaxes(); CL_g_legend(gca(), ["sta1", "sta2"]); // Plot East azimuth <-> elevation for sta1 scf(); for k = 1 : size(tvisi1,2) tk = linspace(tvisi1(1,k), tvisi1(2,k), 100); [az, el] = CL_gm_stationPointing(sta1, CL_interpLagrange(t, pos_ecf, tk)); plot(-az*180/%pi, el*180/%pi); end xtitle("Satellite elevation from sta1 (deg)", "East azimuth (deg)"); CL_g_stdaxes(); | | |
3) Plot ground tracks and visibility "circles"
// ===================================================== // (geocentric) Longitude/latitude plot // ===================================================== scf(); // Plot Earth map CL_plot_earthMap(color_id=color("seagreen")); // Plot ground tracks CL_plot_ephem(pos_ecf, color_id=color("grey50")); // Plot visibility "circles" // min_elev: min elevation for visibility // rmin/rmax: min/max satellite radius (from Earth center) min_elev = 10 * %pi/180; rmin = min(CL_norm(pos_ecf)); rmax = max(CL_norm(pos_ecf)); az = linspace(0,2*%pi,100); CL_plot_ephem(CL_gm_stationVisiLocus(sta1, az, min_elev, rmin), color_id=2); CL_plot_ephem(CL_gm_stationVisiLocus(sta1, az, min_elev, rmax), color_id=2); CL_plot_ephem(CL_gm_stationVisiLocus(sta2, az, min_elev, rmin), color_id=5); CL_plot_ephem(CL_gm_stationVisiLocus(sta2, az, min_elev, rmax), color_id=5); | | |
4) Sun and Moon directions in satellite frame
// ===================================================== // Sun and Moon in Satellite frame // Satellite frame supposed to be "qsw" // ===================================================== M_eci2sat = CL_fr_qswMat(pos_eci, vel_eci); // Satellite->Sun and satellite->Moon directions Sun_dir = M_eci2sat * CL_unitVector(Sun_eci - pos_eci); Moon_dir = M_eci2sat * CL_unitVector(Moon_eci - pos_eci); // Plot angles f=scf(); plot(t-t0, CL_vectAngle(Sun_dir, [0;0;-1])*180/%pi, "r"); plot(t-t0, CL_vectAngle(Moon_dir, [0;0;1])*180/%pi, "b"); xtitle("Angle with satellite frame axis (deg)", "Days"); h = CL_g_select(gca(), "Polyline"); h.thickness = 2; CL_g_legend(gca(), ["Sun <-> -Z", "Moon <-> Z"]); CL_g_stdaxes(); | | |
5) Eclipses
// ===================================================== // Eclipse periods of Sun by Earth // ===================================================== // Radius of Sun and Earth radiusSun = CL_dataGet("body.Sun.eqRad"); radiusEarth = CL_dataGet("body.Earth.eqRad"); // Eclipse indicator (1 = satellite in Earth shadow) ecl = CL_gm_eclipseCheck(pos_eci, Sun_eci, [0;0;0], radiusSun, radiusEarth); // Plot scf(); plot(t-t0, ecl); xtitle("Eclipse periods (1 = satellite in Earth shadow)", "Days"); h = CL_g_select(gce(), "Polyline"); h.thickness = 2; CL_g_stdaxes(); | | |